c-------------------- Linear growth rates  - by V.R.Eke -------------

      function lingro(a,omega,lambda)
c
c     To calculate the linear growth factor D(a) for different cosmological 
c     models. Normalised such that D(1) = 1. (Doesn't include closed models 
c     or lambda models where omega+lambda isn't one.)
c
      implicit none
      real*8 omega,lambda
      real a,func0,func1,func2,x,int,aofx,aofxn
      real xn,w,dn,lingro,sum
      external func0,func1,func2

      w = omega**(-1.0) - 1.0
      sum = omega + lambda
      if (sum.gt.1 .or. omega.le.0 .or.sum.ne.1.and.lambda.gt.0) then
         write(*,*) 'Cannot cope with this cosmology!'
         stop
      endif
      if (omega .eq. 1) then
         lingro = a
      else if (lambda .gt. 0) then
         xn = (2.0*w)**(1.0/3)
         call simp(func2,0.0,xn,int)
         aofxn = ((xn**3.0+2.0)**0.5)*(int/xn**1.5)
         if (a .eq. 1.64) write(*,*) xn,aofxn
         x = a*xn
         call simp(func2,0.0,x,int)
         aofx = ((x**3+2)**0.5)*(int/x**1.5)
         lingro = aofx/aofxn
         if (a .eq. 1.64) write(*,*) x,aofx,lingro
      else
         dn = func1(w)
         x = w*a
         lingro = func1(x)/dn
      endif
 
      end

      function func0(x)
      real x,func0
      func0 = 3/x + (3*((1+x)**0.5)/x**1.5)*log((1+x)**0.5-x**0.5)
      end

      function func1(x)
      real x,func1
      func1 = 1 + 3/x + (3*((1+x)**0.5)/x**1.5)*log((1+x)**0.5-x**0.5)
      end

      function func2(x)
      real x,func2
      func2 = (x/(x**3+2))**1.5
      end

      subroutine simp(func,a,b,s)
      implicit none
      real func,a,b,eps,ost,os,s,st
      integer jmax,j
      external func
      ost = -1.e-30
      os = -1.e-30
      eps = 1.e-5
      jmax = 20
      do j = 1,jmax
         call trapez(func,a,b,st,j)
         s = (4.*st - ost)/3.
         if (abs(s-os) .lt. eps*abs(os)) goto 3
         os = s
         ost = st
      enddo
      pause 'Too many steps.'
3     end

c
c     Trapezium rule from numerical recipes.
c

      subroutine trapez(func,a,b,s,n)
      implicit none
      real func,a,b,s,del,sum,x
      integer j,n,it,tnm
      external func
      save it
      if (n .eq. 1) then
         s = 0.5*(b-a)*(func(a)+func(b))
         it = 1
      else
         tnm = it
         del = (b-a)/tnm
         x = a + 0.5*del
         sum = 0.0
         do j = 1,it
            sum = sum + func(x)
            x = x + del
         enddo
         s = 0.5*(s+(b-a)*sum/tnm)
         it = 2*it
      endif
      end

